
# change directory to the replication location
setwd("D:/Dropbox/jesarey_documents/PS Article on Reviewing")

###########
# combined plot
# unanimity, majority, majority w/ editor, and unilateral editor rules
# 10% reviewer/editor acceptance rate

# 3 reviewers, unanimity required

rm(list=ls())
set.seed(123456)

library(copula)

rho <- seq(from=0, to=0.98, by=0.02)
pr.accept <- 0.1

accept<-c()
for(k in 1:length(rho)){
  reviews <- rCopula(2000, normalCopula(param=c(rho[k], rho[k], rho[k]), dim=3, dispstr="un")) >= (1-pr.accept)
  decisions <- apply(X=reviews, MARGIN=1, FUN=min)
  
  # acceptance rate
  accept[k] <- sum(decisions)/length(decisions)
}

png(file="acceptance-rates.png", width=5, height=6, units="in", res=144)
plot(accept ~ rho, ylim=c(0, 0.13), col=gray(0.5), ylab = "proportion of accepted manuscripts", xlab = "correlation of reviewer opinions", las=2)
y.fit <- predict(loess(accept~rho))
lines(y.fit~rho, lty=3, lwd=2)



# 3 reviewers, majority voting

rm(list=ls())
set.seed(123456)

rho <- seq(from=0, to=0.98, by=0.02)
pr.accept <- 0.1

accept<-c()
for(k in 1:length(rho)){
  reviews <- rCopula(2000, normalCopula(param=c(rho[k], rho[k], rho[k]), dim=3, dispstr="un")) >= (1-pr.accept)
  decisions <- ifelse( apply(X=reviews, MARGIN=1, FUN=sum) >= 2, 1, 0 )
  
  # acceptance rate
  accept[k] <- sum(decisions)/length(decisions)
}

points(accept ~ rho, pch=5, col=gray(0.5))
y.fit <- predict(loess(accept~rho))
lines(y.fit~rho, lwd=2, col=gray(0.5))


# editor as fourth reviewer, 3 out of 4 positive reviews required

rm(list=ls())
set.seed(123456)

rho <- seq(from=0, to=0.98, by=0.02)
pr.accept <- 0.1

accept<-c()
for(k in 1:length(rho)){
  reviews <- rCopula(2000, normalCopula(param=c(rho[k], rho[k], rho[k], rho[k], rho[k], rho[k]), dim=4, dispstr="un")) >= (1-pr.accept)
  decisions.t <- apply(X=reviews, MARGIN=1, FUN=sum)
  decisions <- decisions.t>=3
  
  # acceptance rate
  accept[k] <- sum(decisions)/length(decisions)
}

points(accept ~ rho, pch=4, col=gray(0.5))
y.fit <- predict(loess(accept~rho))
lines(y.fit~rho, lty=2, lwd=2)


# strong editor system with truthful reports

rm(list=ls())
set.seed(123456)

rho <- seq(from=0, to=0.98, by=0.02)
pr.accept <- 0.1

accept<-c()
for(k in 1:length(rho)){
  reviews <- rCopula(2000, normalCopula(param=c(rho[k], rho[k], rho[k], rho[k], rho[k], rho[k]), dim=4, dispstr="un"))
  decisions.t <- apply(X=reviews, MARGIN=1, FUN=mean)
  decisions <- decisions.t >= (1-pr.accept)
  
  # acceptance rate
  accept[k] <- sum(decisions)/length(decisions)
}

points(accept ~ rho, pch=2, col=gray(0.5))
y.fit <- predict(loess(accept~rho))
lines(y.fit~rho, lwd=1)

abline(h=0, col="lightgray")

source(file="LEGEND.r")

# source: http://stackoverflow.com/questions/19053440/r-legend-with-points-and-lines-being-different-colors-for-the-same-legend-item
LEGEND("topleft", bty="n", lty=c(1,2,1,3), pch=c(2,4,5,1), pt.lwd=1, cex=0.8,  lwd=c(1,2,3,2), legend=c("unilateral editor with reviewer advice          ", "majority rule, reviewers + editor          ", "majority rule, reviewers only          ", "unanimity rule, reviewers only          "), col="black", line.col=c("black", "black", gray(0.5), "black"))

dev.off()












# 
# Very large reviewer pool
##########################

rm(list=ls())
library(copula)
set.seed(123456)
rho <- 0.5

accept<-c()

reviews.all <- rCopula(50000, normalCopula(param=rep(rho, 124750), dim=500, dispstr="un"))
reviews <- reviews.all[,1:4]

# strong editor system
decisions.t <- apply(X=reviews, MARGIN=1, FUN=mean)
decisions <- decisions.t >= (1-0.2)
sum(decisions)/length(decisions)

# unanimity voting
decisions.t.2 <- reviews[,1:3] >= (1-0.3)
decisions.2 <- apply(X=decisions.t.2, MARGIN=1, FUN=min)
sum(decisions.2)/length(decisions.2)

# majority voting with active editor
decisions.t.3 <- reviews >= (1-0.2)
decisions.3 <- apply(X=decisions.t.3, MARGIN=1, FUN=sum) >=3
sum(decisions.3)/length(decisions.3)

# majority voting
decisions.t.4 <- reviews[,1:3] >= (1-0.15)
decisions.4 <- apply(X=decisions.t.4, MARGIN=1, FUN=sum) >=2
sum(decisions.4)/length(decisions.4)

reviews.avg <- apply(reviews.all, MARGIN=1, FUN=mean)

png(file="opinion-density.png", width=3.6, height=4.3, units="in", res=144)

plot(density(reviews.avg[which(decisions == T)], from=0.2, to=1), ylim=c(0,6), xlab="mean reader quality evaluation (percentile)", main="Panel 2A: Without Desk Rejection", cex.main=0.9, las=2, cex.lab=0.9)
lines(density(reviews.avg[which(decisions.3 == T)], from=0, to=1), lwd=2, lty=2)
lines(density(reviews.avg[which(decisions.4 == T)], from=0, to=1), lwd=3, col=gray(0.5))
lines(density(reviews.avg[which(decisions.2 == T)], from=0, to=1), lwd=2, lty=3)

legend("topleft", bty="n", lty=c(1,2,1,3), cex=0.6, lwd=c(1,2,3,2), legend=c("unilateral editor with reviewer advice          ", "majority rule, reviewers + editor          ", "majority rule, reviewers only          ", "unanimity rule, reviewers only          "), col=c("black", "black", gray(0.5), "black"))
dev.off()

# proportion of majority voting results at less than the 65th percentile
ecdf(reviews.avg[which(decisions.4 == T)])(0.65)
# proportion of majority voting with active editor results at less than the 65th percentile
ecdf(reviews.avg[which(decisions.3 == T)])(0.65)
# proportion of unanimity voting results at less than the 65th percentile
ecdf(reviews.avg[which(decisions.2 == T)])(0.65)
# proportion of strong editor results at less than the 65th percentile
ecdf(reviews.avg[which(decisions == T)])(0.65)

mean(reviews.avg[which(decisions.4==T)])
mean(reviews.avg[which(decisions.3==T)])
mean(reviews.avg[which(decisions.2==T)])
mean(reviews.avg[which(decisions==T)])

# proportion of majority voting results at less than the 90th percentile
ecdf(reviews.avg[which(decisions.4 == T)])(0.90)
# proportion of majority voting with active editor results at less than the 90th percentile
ecdf(reviews.avg[which(decisions.3 == T)])(0.90)
# proportion of unanimity voting results at less than the 90th percentile
ecdf(reviews.avg[which(decisions.2 == T)])(0.90)
# proportion of strong editor results at less than the 90th percentile
ecdf(reviews.avg[which(decisions == T)])(0.90)

png(file="acceptance-prob.png", width=3.6, height=4.3, units="in", res=144)

# make a plot showing the role of "chance" in the review process
avg.quality <- apply(X=reviews.all, MARGIN=1, FUN=mean)
prob.fit <- loess(decisions.4 ~ avg.quality, degree=0, span=0.05, control=loess.control(surface="interpolate", statistics="approximate", trace.hat="approximate"))
avg.qual.fit <- seq(from=0.2, to=1, by = 0.01)
prob.accept <- predict(prob.fit, newdata=data.frame(avg.quality = avg.qual.fit) )
plot(prob.accept ~ avg.qual.fit, lwd=2, col=gray(0.5), type="l", ylim=c(0, 1), xlab = "mean reader quality evaluation (percentile)", ylab = "probability of acceptance", main="Panel 3A: Without Desk Rejection", yaxs="i", cex.main=0.9, las=2, cex.lab=0.9)


prob.fit <- loess(decisions.2 ~ avg.quality, degree=0, span=0.05, control=loess.control(surface="interpolate", statistics="approximate", trace.hat="approximate"))
avg.qual.fit <- seq(from=0.2, to=1, by = 0.01)
prob.accept <- predict(prob.fit, newdata=data.frame(avg.quality = avg.qual.fit) )
lines(prob.accept ~ avg.qual.fit, lty=3, lwd=2)


prob.fit <- loess(decisions.3 ~ avg.quality, degree=0, span=0.05, control=loess.control(surface="interpolate", statistics="approximate", trace.hat="approximate"))
avg.qual.fit <- seq(from=0.2, to=1, by = 0.01)
prob.accept <- predict(prob.fit, newdata=data.frame(avg.quality = avg.qual.fit) )
lines(prob.accept ~ avg.qual.fit, lty=2, lwd=2)


prob.fit <- loess(decisions ~ avg.quality, degree=0, span=0.05, control=loess.control(surface="interpolate", statistics="approximate", trace.hat="approximate"))
avg.qual.fit <- seq(from=0.2, to=1, by = 0.01)
prob.accept <- predict(prob.fit, newdata=data.frame(avg.quality = avg.qual.fit) )
lines(prob.accept ~ avg.qual.fit)

legend("topleft", bty="n", cex=0.6, lty=c(1,2,1,3), lwd=c(1,2,3,2), legend=c("unilateral editor with reviewer advice        ", "majority rule, reviewers + editor        ", "majority rule, reviewers only        ", "unanimity rule, reviewers only         "), col=c("black", "black", gray(0.5), "black"))

dev.off()



# impose desk rejections for the papers the editor thinks are the worst
not.desk <- which(reviews[,4] >= 0.5)

# strong editor system
decisions.t <- apply(X=reviews[not.desk,], MARGIN=1, FUN=mean)
decisions <- decisions.t >= (1-0.13)
sum(decisions)/length(decisions)

# unanimity voting
decisions.t.2 <- reviews[not.desk,1:3] >= (1-0.21)
decisions.2 <- apply(X=decisions.t.2, MARGIN=1, FUN=min)
sum(decisions.2)/length(decisions.2)

# majority voting with active editor
decisions.t.3 <- reviews[not.desk,] >= (1-0.15)
decisions.3 <- apply(X=decisions.t.3, MARGIN=1, FUN=sum) >=3
sum(decisions.3)/length(decisions.3)

# majority voting
decisions.t.4 <- reviews[not.desk,1:3] >= (1-0.09)
decisions.4 <- apply(X=decisions.t.4, MARGIN=1, FUN=sum) >=2
sum(decisions.4)/length(decisions.4)

reviews.avg <- apply(reviews.all[not.desk,], MARGIN=1, FUN=mean)

png(file="opinion-density-desk.png", width=3.6, height=4.3, units="in", res=144)

plot(density(reviews.avg[which(decisions == T)], from=0.2, to=1), ylim=c(0,6), xlab="mean reader quality evaluation (percentile)", main="Panel 2B: With Desk Rejection", cex.main=0.9, las=2, cex.lab=0.9)
lines(density(reviews.avg[which(decisions.3 == T)], from=0, to=1), lwd=2, lty=2)
lines(density(reviews.avg[which(decisions.4 == T)], from=0, to=1), lwd=3, col=gray(0.5))
lines(density(reviews.avg[which(decisions.2 == T)], from=0, to=1), lwd=2, lty=3)

legend("topleft", bty="n", cex=0.6, lty=c(1,2,1,3), lwd=c(1,2,3,2), legend=c("unilateral editor with reviewer advice          ", "majority rule, reviewers + editor          ", "majority rule, reviewers only          ", "unanimity rule, reviewers only          "), col=c("black", "black", gray(0.5), "black"))

dev.off()

# proportion of majority voting results at less than the 65th percentile
ecdf(reviews.avg[which(decisions.4 == T)])(0.65)
# proportion of majority voting with active editor results at less than the 65th percentile
ecdf(reviews.avg[which(decisions.3 == T)])(0.65)
# proportion of unanimity voting results at less than the 65th percentile
ecdf(reviews.avg[which(decisions.2 == T)])(0.65)
# proportion of strong editor results at less than the 65th percentile
ecdf(reviews.avg[which(decisions == T)])(0.65)

mean(reviews.avg[which(decisions.4==T)])
mean(reviews.avg[which(decisions.3==T)])
mean(reviews.avg[which(decisions.2==T)])
mean(reviews.avg[which(decisions==T)])

# proportion of majority voting results at less than the 90th percentile
ecdf(reviews.avg[which(decisions.4 == T)])(0.90)
# proportion of majority voting with active editor results at less than the 90th percentile
ecdf(reviews.avg[which(decisions.3 == T)])(0.90)
# proportion of unanimity voting results at less than the 90th percentile
ecdf(reviews.avg[which(decisions.2 == T)])(0.90)
# proportion of strong editor results at less than the 90th percentile
ecdf(reviews.avg[which(decisions == T)])(0.90)

png(file="acceptance-prob-desk.png", width=3.6, height=4.3, units="in", res=144)

# make a plot showing the role of "chance" in the review process
avg.quality <- reviews.avg
prob.fit <- loess(decisions.4 ~ avg.quality, degree=0, span=0.05, control=loess.control(surface="interpolate", statistics="approximate", trace.hat="approximate"))
avg.qual.fit <- seq(from=0.2, to=1, by = 0.01)
prob.accept <- predict(prob.fit, newdata=data.frame(avg.quality = avg.qual.fit) )
plot(prob.accept ~ avg.qual.fit, lwd=2, col=gray(0.5), type="l", ylim=c(0, 1), xlab = "mean reader quality evaluation (percentile)", ylab = "probability of acceptance", main="Panel 3B: With Desk Rejection", yaxs="i", cex.main=0.9, las=2, cex.lab=0.9)


prob.fit <- loess(decisions.2 ~ avg.quality, degree=0, span=0.05, control=loess.control(surface="interpolate", statistics="approximate", trace.hat="approximate"))
avg.qual.fit <- seq(from=0.2, to=1, by = 0.01)
prob.accept <- predict(prob.fit, newdata=data.frame(avg.quality = avg.qual.fit) )
lines(prob.accept ~ avg.qual.fit, lty=3, lwd=2)


prob.fit <- loess(decisions.3 ~ avg.quality, degree=0, span=0.05, control=loess.control(surface="interpolate", statistics="approximate", trace.hat="approximate"))
avg.qual.fit <- seq(from=0.2, to=1, by = 0.01)
prob.accept <- predict(prob.fit, newdata=data.frame(avg.quality = avg.qual.fit) )
lines(prob.accept ~ avg.qual.fit, lty=2, lwd=2)


prob.fit <- loess(decisions ~ avg.quality, degree=0, span=0.05, control=loess.control(surface="interpolate", statistics="approximate", trace.hat="approximate"))
avg.qual.fit <- seq(from=0.2, to=1, by = 0.01)
prob.accept <- predict(prob.fit, newdata=data.frame(avg.quality = avg.qual.fit) )
lines(prob.accept ~ avg.qual.fit)

legend("topleft", bty="n", cex=0.6, lty=c(1,2,1,3), lwd=c(1,2,3,2), legend=c("unilateral editor with reviewer advice        ", "majority rule, reviewers + editor        ", "majority rule, reviewers only        ", "unanimity rule, reviewers only         "), col=c("black", "black", gray(0.5), "black"))

dev.off()






###########
# plot average review quality against percentile
# strong editor system with truthful reports
# two subfield discipline

rm(list=ls())
set.seed(123456)

library(copula)

N <- 250              # size of each subdiscipline

rho.hi <- 0.9         # within-subfield correlation
rho.lo <- 0.1         # between-subfield correlation

rho.vec <- c()
for(i in (N-1):0){
  rho.vec <- c(rho.vec, rep(rho.hi, i), rep(rho.lo, N))
}
rho.vec <- c(rho.vec, rep(rho.hi, sum(1:(N-1))))

reviews.all <- rCopula(50000, normalCopula(param=rho.vec, dim=(N*2), dispstr="un"))
# calculate average correlation
cor.mat <- cor(reviews.all)
mean(cor.mat)

reviews <- reviews.all[,c(1:2, 499:500)]

# strong editor system
decisions.t <- apply(X=reviews, MARGIN=1, FUN=mean)
decisions <- decisions.t >= (1-0.22)
sum(decisions)/length(decisions)

review.avg <- apply(X=reviews, MARGIN=1, FUN=mean)
rank.quant <- ecdf(review.avg)(review.avg)
o<-order(rank.quant)

png(file="opinion-density-subdis.png", width=3.6, height=4.3, units="in", res=144)

plot(density(review.avg[which(decisions == T)], from=0.2, to=1), ylim=c(0,13.5), xlab="mean reader quality evaluation (percentile)", main="Panel 4A: Without Desk Rejection", cex.main=0.9, las=2, cex.lab=0.9)

mean(review.avg[which(decisions==T)])


# repeat with one field only, low correlation
rm(list=ls())
set.seed(123456)
rho <- 0.5

accept<-c()

reviews.all <- rCopula(50000, normalCopula(param=rep(rho, 124750), dim=500, dispstr="un"))
reviews <- reviews.all[,1:4]

# strong editor system
decisions.t <- apply(X=reviews, MARGIN=1, FUN=mean)
decisions <- decisions.t >= (1-0.2)
sum(decisions)/length(decisions)

review.avg <- apply(reviews.all, MARGIN=1, FUN=mean)

lines(density(review.avg[which(decisions == T)], from=0.2, to=1), lty=2)

mean(review.avg[which(decisions==T)])




# repeat with one field only, high correlation
rm(list=ls())
set.seed(123456)
rho <- 0.75

accept<-c()

reviews.all <- rCopula(50000, normalCopula(param=rep(rho, 124750), dim=500, dispstr="un"))
reviews <- reviews.all[,1:4]

# strong editor system
decisions.t <- apply(X=reviews, MARGIN=1, FUN=mean)
decisions <- decisions.t >= (1-0.15)
sum(decisions)/length(decisions)

review.avg <- apply(reviews.all, MARGIN=1, FUN=mean)

lines(density(review.avg[which(decisions == T)], from=0.2, to=1), lwd=2)

mean(review.avg[which(decisions==T)])




legend("topleft", bty="n", cex=0.6, lty=c(1,2,1), lwd=c(1,1,2), col=c("black", "black"), legend=c("two subdisciplines, rho = 0.9 or 0.1          ", "all readers correlated at rho = 0.5          ", "all readers correlated at rho = 0.75          ") )


dev.off()




###########
# plot average review quality against percentile
# strong editor system with truthful reports
# two subfield discipline
# with desk rejection

rm(list=ls())
set.seed(123456)

library(copula)

N <- 250              # size of each subdiscipline

rho.hi <- 0.9         # within-subfield correlation
rho.lo <- 0.1         # between-subfield correlation

rho.vec <- c()
for(i in (N-1):0){
  rho.vec <- c(rho.vec, rep(rho.hi, i), rep(rho.lo, N))
}
rho.vec <- c(rho.vec, rep(rho.hi, sum(1:(N-1))))

reviews.all <- rCopula(50000, normalCopula(param=rho.vec, dim=(N*2), dispstr="un"))
reviews <- reviews.all[,c(1:2, 499:500)]

# impose desk rejection
not.desk <- which(reviews[,4] >= 0.5)

# strong editor system
decisions.t <- apply(X=reviews[not.desk,], MARGIN=1, FUN=mean)
decisions <- decisions.t >= (1-0.15)
sum(decisions)/length(decisions)

review.avg <- apply(X=reviews[not.desk,], MARGIN=1, FUN=mean)
rank.quant <- ecdf(review.avg)(review.avg)
o<-order(rank.quant)

png(file="opinion-density-subdis-desk.png", width=3.6, height=4.3, units="in", res=144)

plot(density(review.avg[which(decisions == T)], from=0.2, to=1), ylim=c(0,13.5), xlab="mean reader quality evaluation (percentile)", main="Panel 4B: With Desk Rejection", cex.main=0.9, las=2, cex.lab=0.9)

mean(review.avg[which(decisions==T)])


# repeat with one field only, low correlation
rm(list=ls())
set.seed(123456)
rho <- 0.5

accept<-c()

reviews.all <- rCopula(50000, normalCopula(param=rep(rho, 124750), dim=500, dispstr="un"))
reviews <- reviews.all[,1:4]

# impose desk rejection
not.desk <- which(reviews[,4] >= 0.5)

# strong editor system
decisions.t <- apply(X=reviews[not.desk,], MARGIN=1, FUN=mean)
decisions <- decisions.t >= (1-0.13)
sum(decisions)/length(decisions)

review.avg <- apply(reviews.all[not.desk,], MARGIN=1, FUN=mean)

lines(density(review.avg[which(decisions == T)], from=0.2, to=1), lty=2)

mean(review.avg[which(decisions==T)])


# repeat with one field only, higher correlation
rm(list=ls())
set.seed(123456)
rho <- 0.75

accept<-c()

reviews.all <- rCopula(50000, normalCopula(param=rep(rho, 124750), dim=500, dispstr="un"))
reviews <- reviews.all[,1:4]

# impose desk rejection
not.desk <- which(reviews[,4] >= 0.5)

# strong editor system
decisions.t <- apply(X=reviews[not.desk,], MARGIN=1, FUN=mean)
decisions <- decisions.t >= (1-0.09)
sum(decisions)/length(decisions)

review.avg <- apply(reviews.all[not.desk,], MARGIN=1, FUN=mean)

lines(density(review.avg[which(decisions == T)], from=0.2, to=1), lwd=2)

mean(review.avg[which(decisions==T)])



legend("topleft", bty="n", cex=0.6, lty=c(1,2,1), lwd=c(1,1,2), col=c("black", "black"), legend=c("two subdisciplines, rho = 0.9 or 0.1          ", "all readers correlated at rho = 0.5          ", "all readers correlated at rho = 0.75          ") )

dev.off()





